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Zigo Haras and Shlomo Ta’asan * 

Department of Applied Mathematics and Computer Science 
The Weizmann Institute of Science 
and 

Institute for Computer Applications in Science and Engineering 


Abstract 

Finite difference schemes for the evaluation of first and second derivatives are presented. 
These second order compact schemes were designed for long-time integration of evolution equa- 
tions by solving a quadratic constrained minimization problem. The quadratic cost function 
measures the global truncation error while taking into account the initial data. The resulting 
schemes are applicable for integration times fourfold, or more, longer than similar previously 
studied schemes. A similar approach was used to obtain improved integration schemes. 
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1 Introduction 


The simulation of hyperbolic partial differential equations often requires long-time integration. 
The physical phenomena described by these equations typically possess a range of space and time 
scales, turbulent fluid flow is a common example. Accurate numerical simulation of this type of 
processes requires proper representation of all the relevant physical scales in the numerical model. 
These requirements lead recently to new interest in Pade approximations also known as compact 
finite difference schemes [7]. 

Compact finite difference schemes had long been known and used in numerical analysis [1, 2, 3]. 
They offer a means of obtaining high order approximations to differential operators using narrow 
stencils. This is achieved by treating the sought derivatives as unknowns and solving a system of 
equations for them. Typically, the resulting matrices are tridiagonal or pentadiagonal, which can 
be efficiently solved . 

In [7] a class of highly accurate compact schemes for first, second and higher derivatives were 
presented and analyzed. A notion of resolving efficiency was introduced which should measure 
the accuracy with which the finite difference approximation represents the exact solution over 
the full range of length scales that can be realized on a given grid. This criterion was then used 
to compare various schemes and motivated the design of a new class of schemes, the so-called 
schemes with spectral like resolution . These are fourth order pentadiagonal systems with seven 
points stencils. Their improved resolution characteristics were obtained by giving up on high 
formal accuracy; instead, requiring that the symbol of the discrete difference operator should 
agree with the differential operator at three prescribed high frequencies. However, the resolving 
efficiency is a too crude measure as it assumes that all frequencies occur with similar magnitude in 
the initial data. In the present paper it is shown that for problems with various initial conditions 
these schemes are far from optimal. 

A different measure for evaluating finite difference schemes is the L 2 norm of the local trunca- 
tion error. This measure which takes into account the Fourier components present in the solution 
and their amplitude was applied in [9] to design explicit time marching schemes (i.e., discretiz- 
ing time and space simultaneously) by solving analytically constrained minimization problems 
with quadratic cost. This error measure seems more adequate for comparing difference schemes. 
However, the simultaneous treatment of time and space results in very complex optimization 
problems. The generalization of the approach of [9] to problems in higher dimensions requires 
solving nonlinear constrained optimization over a large set of parameters. It is also hard to apply 
that approach to compact schemes. This complexity and the use of analytic rather than numerical 
methods makes the suggested approach impractical. 

In [4] a heuristic derivation was done by minimizing the weighted error (in the Fourier space) 
of the discrete and continuous operators. 

The present paper uses the same cost function as [9] with several important differences. First, 
improved bounds on the truncation error were derived. These enabled us to treat time and 
space discretizations separately. This greatly simplifies the minimization problem by reducing 
it to two lower dimensional problems. Further simplification is obtained by optimizing each 
partial derivative separately, instead of approximating the whole differential operator as was done 
in [9]. These reductions of problem complexity resulted in a simple and general approach to 
synthesis of discretization schemes. It enabled us to design highly accurate compact difference 
schemes and integration formulas for various operators and initial data. The resulting second 
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order approximations proved to be robust to perturbations in the spectrum of the initial data, 
exhibiting resolution superior to other known schemes. 

The organization of the paper is as follows. In Section 2 Fourier analysis is used to obtain 
bounds on the truncation error. In section 3 approximations to derivatives are presented, for the 
first and second derivatives and first derivative at mid cell points. Appendices A-C list coefficients 
for these derivatives for various stencils and initial conditions. Improved time integration schemes 
are developed in section 4, and their coefficients are listed at Appendix D. Section 5 discusses 
generalization of the present approach to more complex problems. Numerical results are presented 
in section 6. Concluding remarks are made in section 7. 


2 Bounds on the truncation error 


The application of Fourier analysis for the design and evaluation of finite difference schemes can 
be found in many sources, e.g., [9, 10). In [12] the use of Fourier analysis in the numerical 
approximation of hyperbolic problems is extensively discussed. 

In the following section bounds on the L 2 norm of the error in the discrete solution are derived, 

accounting for the effect of discretization both in space and time. These estimates are used in 

subsequent sections to design highly accurate schemes. 

Consider a linear constant coefficient partial differential equation with periodic boundary 
conditions of the form : 

% = < 2 -’> 

u(x,0) = uq(x) ( 2 - 2 ) 


Further assume that the solution of equation (2.2) does not grow in time. 

The discrete analog of this equation can be written as : 

(2.3) 

(2.4) 


< +1 = P(fc,A*K 

u° h = Uq 


where h is the meshsize in space, and P(h, At) is a stable finite difference approximation. 

We would like to bound the L 2 norm of the error in the discrete solution, for the initial value 
tio, given by : 

e 2 («Af ; u 0 ) = ||u(nAt) - <|| 2 = ||e LnA ‘«o - P n (h, A*)«o|| 2 (2.5) 


In the sequel t will be used instead of nAt to simplify notation. 

(2.5) yields 

J ^ - P^wfe)! |«o(u>/i)| 2 du 


The Fourier transform of eq. 
< ( 2 - 6 ) 


/* _ e L h (»h)t | + |e ih ( wA )‘ - P n (u>h) \f |tio(wh)| 2 <iw 


(2.7) 


Where L h is the discrete operator approximating L, and L, L h are their corresponding symbols. 
Thus, the space and time discretizations errors can be bounded separately. 
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Denote L(uk) = L R {uh ) + L^uh ) for the real and imaginary parts of L(uh), respectively; 


and use a similar decomposition for L h (u)h). Then 

e L(wh)t _ e L h (wh)t _ e L(wh)t ^ _ e (L*(ujh)-Li(wh))t e (L%(u,h)-L R (wh.))t'j (2.8) 

The assumption that the solution does not grow in time implies that 

| e ^M)‘| <1 t > 0 (2.9) 

Assume that the discrete solution does not grow in time, either. That is 

| e L h (^)i|< 1 *> 0 (2.10) 

For real numbers 0, a with a < 0, simple geometric considerations yield the following bound 

|1 -e'V| < |0| + |l-e“| (2.11) 

Combining bounds (2.9) and (2.11) and assuming L^(u>K) — L R (u>h ) < 0 results in : 

\ e L^h)t _ e L»(»K)t\ < _i / ( w fc)|t + |l-e^" fc >- X ’*< wA))t | (2.12) 

< |L/(w/i) — Lj(u?h)\t + \L R (u>h) — L R (u>h)\t (2.13) 


If ~L h R {uh) - L R (u>h) > 0, this bound can be obtained using the same argument when factoring 
e L h (wh)t in (2. 8 ). 

Denote by e(t; uo) the error due to spatial discretization only, when the initial data is uq. For 
a final time T, using (2.13) yields : 

e 2 (T; «o) <T 2 (| L)(uh) - ti(uh)\ + | L h R (uh) - L R (ujh ) |) 2 \u 0 (ujh)\ 2 du (2.14) 

7T 

Therefore, a difference scheme minimizing the integral (2.14) with respect to initial value uo 
will better resolve, in the L 2 norm sense, the frequencies occurring in the solution. 

The time integration operator satisfies : 

e L h (wh)nAt _ pn( w /j) — f e L h (wK)&t _ p^hf) Y") g£ h (^)jAt pn-l-j^^) (2.15) 

J=1 

Under the previous assumptions 

| e L h (w/i)Ai| < x ( 2 . 16 ) 

\P(e>h)\ < 1 ( 2 - 17 ) 

Therefore 

| ^ e i' l (u/h)iAtpn-l- J '( u; ^)j < Cn (2.18) 

i = i 

where C = 0(1). Hence 

| e Lh W - P n (uh)\ < _ p(uh ) | (2.19) 

v 
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Denote by e(t ; uq) the error due to time discretization only, when the initial value is uq. For a 
final time T, the bound (2.19) implies : 


e 2 (T,« 0 ) < (|j)7 * " n^)i 2 |«oM)| 2 du; 
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( 2 . 20 ) 


Combining of these estimates yields a bound on the L 2 norm of the global truncation error : 


e 2 (T; tie) < T 2 J\ [| t)(uh) - L t (uh)\ + | L%(u>h) - i R (uh)\ + £|e^ (wfc)At - P(uh)\ 


\uo{uh)\ 2 du 
( 2 . 21 ) 


These estimates will be used in the following sections to design improved compact schemes and 
integration formulas. 


3 Approximating derivatives 


3.1 Approximation of the first derivative 

Consider a uniformly spaced mesh whose nodes are indexed by i and its meshsize is given by 
h = jj, where N + 1 is the number of grid points. The variable at node i is £,• = ih and the 
function value at the nodes, /, = /(x,), are given for 0 < i < N . An approximation f' to the first 
derivative ^ (x,) should be computed as a linear combination of the function values at neighboring 
grid points. Compact finite difference schemes regard the approximation J\ as unknown and a 
system of equations is solved to approximate the first derivative at all nodes, simultaneously. 
Thus, unlike in finite differences, the derivative at node i depends on function values at all other 
nodes. 

Following [7] we use approximations of the the form : 


Pfl - 2 + “//- 1 + f'i + a fi + 1 + Pfl+2 ~ c 


fi + 3 - ft— 3 . ,fi+ 2 - ft— 2 . _ fi+1 - ft - 1 


6h 


+ b 


4h 


+ a- 


2h 


(3.1) 


A second order approximation can be obtained by adding a constraint that the Taylor expansion 
on both sides should agree until the second order term. The following relation must hold : 


a + 6 + c — l - ! - 2ft T 2 /? 


(3.2) 


Higher order schemes may be obtained by further matching the next terms in the expansion [7]. 
However, in this paper merely second order accuracy is enforced. 

The symbol of the differentiation operator is given by : 

L(uh) — iuh (3-3) 


Whereas the symbol of the discrete approximation (3.1) is 

i k (uh) = i 


.a sm(u>h) + | sin(2u>h) + | sin(3u>/i) 


1 + 2ft cos(o»/i) -f 2/? cos(2u>h) 


(3.4) 
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In view of the bound (2.14), define the following constrained minimization problem which its 
solution should yield a compact scheme with improved resolution properties : 

min [ h \t h (ujh) — L(uh)\ \u(uh)\ 2 du (3*5) 

a,6,c,ar t /? I * 

under the constraint 

a+b+c= l + 2a + 2/? (3.6) 

where L(u>h) and L k (uh) are given by (3.3) and (3.4), respectively. 

Although the problem was formulated as a constrained minimization problem, it can be trans- 
formed by substitution to an unconstrained minimization problem on a reduced set of parameters. 
Moreover, setting some of the parameters to zero further reduces the dimension of the problem. 
Since tridiagonal systems of equations are more amenable to numerical solution than pentadi- 
agonal ones, setting (3 = 0 seems a plausible choice. Similar considerations might suggest using 
a narrower stencil obtained by setting c = 0, as well. All those possibilities are presented in 
Section 6, and several sets of coefficients for different initial data are listed in Appendix A. 


3.2 Approximation of the second derivative 

The derivation of compact schemes for the second derivative proceeds in an analogous way to the 
first derivative. The starting point is an approximation of the form 


+ 4 '+< + , + MU = + » / * a ~w + / - + 

where f- 1 is the approximation to the second derivative at node i. Matching the Taylor series 
coefficients on both sides of (3.7) yields condition (3.2) for the second order accuracy. 

The symbol of the second derivative is given by : 


L(ujh) = —u> 2 h 2 


(3.8) 


The symbol of the discrete approximation (3.7) is : 


l h {uh) = - 


2 a (cos(a;/t) — 1) + f (cos(2a >h) — 1) + (cos(3 uh) — 1) 

1 + 2a cos(wh) + 2/3 cos(2wh) 


(3.9) 


The constrained minimization problem which solution is the sought scheme can be formulated 

as : w 

min /* \L\uh) - L(wh)| 2 \u(ujh)\ 2 du> (3.10) 

a,b,c,a,(} J-lfc I I 

under the constraint 

a + b + c=l + 2a + 2(3 (3.11) 

Now, however, L(u>h ) and L h (uh) are given by (3.8) and (3.9), respectively. 


5 



3.3 Approximating first derivative on a cell centered mesh 

The approximation of the first derivative at the cell centered mesh is of the form : 


/,-+5 - fi-k fi+k - fi-1 Ji + 

flf'i - 2 + a f'-i + fi + a /.'+i + Pfi+2 = c ' - + b ; h a 




5/i 


3ft 


ft 


(3.12) 


The second order of the approximation is guaranteed by condition (3.2). 
The symbol of the differentiation operator is : 


L(u>h) — iuh 


(3.13) 


While the symbol of the discrete approximation (3.12) is : 


.. , 2.sm(f)+f sm(M) + y sin (5a>) 

^ 1 1 + 2 a cos(c oh) + 2/3 cos(2u )h) 


(3.14) 


A constrained minimization problem of the same type as in the previous sections was formulated 
and solved for these symbols. 


4 Approximation of the integration operator 

The design of integration schemes is substantially limited by the stability requirement which 
renders high order schemes computationally costly. It is well known [5] that an explicit k ih order 
Runge-Kutta method, with k > 4, should employ at least fc+1 stages or function evaluations . This 
large number of function evaluations makes high order schemes impractical. Therefore, efforts 
have been made to obtain schemes of lower order with improved characteristics. Within this 
approach, the free variables in the Runge-Kutta schemes were set to yield better truncation error 
[5] or extended stability region [6]. A generalization of this idea is to give up on formal accuracy 
in order to obtain better approximation of the wavenumbers relevant to the problem solved. 

The discrete time integration of linear constant coefficient partial differential equation 

d i - < 41 > 

amounts to approximation of the exact discrete solution e Lht Uo. Therefore, the integration scheme 
may be written as 

p n (L h At) = j2 a ^ LhAt y ( 4 - 2 ) 

t=0 

where a* may depend on L h . The order of the integration scheme is determined by the number 
of first terms a; which agrees with the Taylor expansion of e x . 

The derivation of the integration schemes is similar to that of derivative discretization, i.e,, 
a constrained quadratic optimization problem is formulated based on the error estimate (2.20). 
The solution of this minimization problem yields an improved integration scheme. However, the 
derivation of integration schemes is more involved than the generation of compact schemes since 
the stability condition leads to a nonlinearly constrained minimization problem. 
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Following (2.20) the next optimization problem is defined : 

min \e Lh ^ Ai - P n {L h {uh)At)\ 2 \u{uh)\ 2 du (4.3) 

°* J ~l 

subject to the constraints 

ai = ^ 0<t<p (4.4) 

l\ 

\P n (L k {uh)At )\ 2 < 1 (4-5) 

where L h is the discrete approximation of L and p is the order of the n stage formula. Condition 
(4.4) can be treated by substitution, but the stability condition requires an explicit treatment. 

In accordance with our general approach, we believe that second order formal accuracy suffices. 
It remains to determine the number of stages in the integration formula. This should be chosen 
to assure that the error in space and time discretizations (2.14) and (2,20), respectively, will be 
of similar magnitude. In the present work five stage schemes of second order were investigated, 
i.e., n = 5 and p = 2. Integration formulas were obtained for optimized seven point tridiagonal 
compact schemes approximating the first derivative, and were tested for the advection equation 
in one and two space dimensions. 

An important feature of the present approach is that once a feasible minimum has been found 
for a prescribed initial value and a given CFL number, the generated scheme will be stable for 
this data. This might enable the use of somewhat larger time steps. 

5 Approximation of differential operators 

The method introduced in the previous sections for generating optimal finite difference approxi- 
mations for derivatives and time integration schemes for linear constant coefficient equations can 
be extended to more general equations. These ideas can be easily adapted to handle with similar 
efficiency more involved problems. 

The error bounds derived in Section 2 can be generalized for d dimensional problems; noting 
that the same proof holds for the d dimensional case after changing the integration over [— f , f ] 
to multi integration over the box ] d . This suggests that approximation of the differen- 

tial equation should be obtained by solving constrained optimization problems in d dimensional 
Fourier space for a large set of parameters. For some equations, solving this large minimization 
problem might be essential to achieve accurate schemes. Quite often, though, a set of simpler 
minimization problems can be obtained by optimizing each partial derivative separately, resulting 
in highly accurate approximations. 

The approach, which was successfully tested in the present paper, divides the optimization 
process into two stages. In the first stage a set of schemes are designed for a large enough 
variety of typical initial data (e.g., Gaussians with different parameters, in our examples). This 
precomputation is performed once and its results are used in subsequent simulations. In the 
second stage, the actual simulation, the initial data uq is Fourier transformed to obtain uo. 
The discretization of the partial derivatives is determined by approximating fio as a product of 
one dimensional functions for which optimized schemes were designed. Each partial derivative 
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is discretized using the corresponding one dimensional optimized scheme. The time marching 
scheme is selected from the set of scheme corresponding to the approximating one dimensional 
functions. In the present work the selection was done by computing the L 2 error norm of each 
candidate integration scheme when applied to the approximate initial data with the already 
determined discretizations, and selecting the minimum norm scheme. This computation, too, can 
be done prior to the actual simulation for a large set of approximated initial data. Thus, the 
marching scheme selection can be done by looking up in a precomputed table. The robustness of 
the proposed schemes to perturbations in the initial data yields this optimization very efficient; 
as can be seen in the numerical results presented in Section 6. It should be noted that the 
time required to obtain an appropriate scheme using this approach is negligible relative to the 
simulation time. 

When the frequencies present in the solution change with time, e.g., due to time dependent 
source term, the computation of the optimized schemes should be repeated once a large cumulative 
change has occurred. Still, -the relative cost of of this computation is minimal. 

The Fourier transform gives the energy content of the whole initial data. It may occur, that 
the initial data is smooth at some regions of the computational domain and oscillatory in others, 
in which case the designed approximation will give good performances over the whole domain. 
One can do better by computing a different scheme for each region and using a weighted sum 
of the resulting schemes near region boundaries. This requires computing the Fourier transform 
locally in each region. The localization to a particular region can be achieve by multiplying uq by 
a C°° function with a compact support which encloses the region. 

In some cases, systems of equations may be treated in a similar way. Look first at a one 
dimensional first order system 

u t = Au x (5.6) 

u(x,0) = uo(x) 

where A is a p x p symmetric matrix. Let A = P~ l AP be a diagonal matrix, and denote v = Pu . 
The discretization of the system 

v t = Av x (5.7) 

t;(x,0) = Puq(x) 

can be done in an analogous way to the scalar case, except for the time marching scheme which 
should be chosen from a set of candidate schemes (as for the multidimensional scalar equations). 
Thus, highly accurate discretization of the system (5.6) can be achieved by first discretizing (5.7) 
and using the identity u x — P~ l v x . For systems in higher dimension 

u < = ( 5 - 8 ) 

i=l ° X% 

u(z,0) = u 0 (x) 

each partial derivative should be optimized separately. In this case, we require that all Ai will be 
symmetric, but it is not necessary that they are simultaneously diagonalizable. 

The proposed schemes might be useful for nonlinear equations, as well. There, one should 
design the schemes for the linearized equation; and will be obliged modify them once a large 
change in the amplitude of the wavenumbers appearing in the solution occur. 
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6 Numerical results 

6*1 Approximation of derivatives 

The constrained minimization problem (3.5) for the space discretization can be easily solved by 
substitution using (3.2). Differentiation of the resulting quadratic form provides a set of necessary 
conditions holding at the minimum. This nonlinear system can be solved using Newton method, 
yielding a local minimum. Since the schemes obtained using this process significantly improve 
previously known schemes [7], no attempts were made to find the other zeroes of the nonlinear 
system, searching for better minima. 

Three types of schemes were studied : (a) tridiagonal with five points stencil, i.e., (3 = c = 0 , 
(b) tridiagonal with seven points stencil, i.e., (3 = 0, (c) pentadiagonal with seven points stencil. 
The initial approximation to the Newton iteration was, typically, a compact scheme with the same 
structure, taken from [7]. 

It can be observed, in figures 1 and 5, that the modulos of symbol of the optimized penta- 
diagonal scheme for the first and second derivative is larger then the modulus of the differential 
symbol. This error is exceedingly larger for schemes generated to approximate narrower spectra. 
The overshooting occurs in the highest end of the spectrum for wavenumbers not appearing in 
the solution. However, since the stability of a scheme is determined by the values assumed by 
L h (uh) [11], this type of scheme is applicable only with small CFL. Moreover, the desired ro- 
bustness is limited by this phenomenon. Therefore, this behavior of the approximation can not 
be ignored. A possible remedy can be found in searching for other minimizers of the quadratic 
form. Using the tridiagonal scheme as initial approximation for the Newton process converged 
to solutions without this limiting property but with reduced resolution, similar to the tridiago- 
nal schemes. Other possible directions, e.g., further looking for other minima or penalizing in 
the cost function for this behavior were not explored. This is since we believe that for practical 
applications pentadiagonal systems are too costly to solve, whereas the tridiagonal schemes offer 
similar resolution characteristics, are easier to solve and do not suffer from this deficiency. The 
pentadiagonal scheme are given mainly for theoretical reasons as a counterpart to the spectral-like 
approximations. 

A proper appreciation of the superiority of the proposed schemes can be gained by integrating 
with them hyperbolic equations for long time, provided the integration process introduces only 
negligible numerical errors. This requirement necessitates either using high order integration 
schemes or employing exact integration, as was done in the present work. The experiments 
described in next subsections clearly demonstrate the superior behavior of the proposed optimized 
schemes. 

6.1.1 Approximation of the first derivative 

2 

Compact finite difference schemes were designed and tested for initial data of the form e~ au; for 
several values of a. In Figure 1, the symbols of schemes corresponding to a = 2 are plotted, as 
well as the weighted error 

\L((jh) — L h (u>h)\\u(u>h)\ (6.9) 
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for the more accurate schemes. The coefficients of the optimized schemes can be found in Appendix 
A. The coefficients of the other schemes were taken from [7]. For scheme (a) the coefficients were 

a = l /3 = 0,a = y,6=^c = 0 (6.10) 

The coefficients of scheme (c) were 

a = -,/S = 0,a=-,6=-,c=-- (6.11) 

The coefficients of the spectral-like scheme (e) were 

a = 0.5771439, (3 = 0.0896406, a = 1.3025166,6 = 0.99355, c = 0.03750245 (6.12) 


It can be seen that each optimized scheme better approximates the differential operator than 
its non-optimized counterpart. In Figure 1, one can observe that although the symbol of the 
spectral-like pentadiagonal scheme follows the differential symbol for more wavenumbers than the 
tridiagonal scheme, the L 2 norm of truncation error of tridiagonal scheme is somewhat smaller 
for this data. This can be explained by noting that the error of the tridiagonal scheme is mainly 
in the high frequencies while the spectral-like scheme has large error at the smoother Fourier 
components where the present initial data has more energy. The spectral-like scheme attains 
better resolution at the expense of larger error in lower frequencies. The error in the optimized 
schemes is significantly smaller than in their counterparts. More precisely, computing the error 
norms reveals that the error in the tridiagonal scheme is about six times larger than in the 
optimized tridiagonal scheme while error norm of the spectral-like scheme is about seventeen 
times larger than in the optimized pentadiagonal. The plot of the absolute value of the error 
also reveals that the L 2 norm was used as a minimization criteria. This can be seen from the 
several sign changes of the error of the optimized schemes, being in accordance with the averaging 
property of the chosen norm. 

Figure 2 demonstrates the better resolution of the optimized scheme by exact integrating in 
Fourier space on a 32 points grid with the pentadiagonal spectral-like scheme and the pentadiag- 
onal optimized scheme the equation 


a = 0.8 (a being the CFL number) was used. It is shown that at time T = 10000, the error in the 
solution using the optimized scheme is smaller than the error at time T = 1000 when using the 
spectral-like scheme. This suggests that the optimized scheme can be used for integration time 
at least ten times longer than the spectral-like scheme, in close accordance with the ratio of the 
error norms. 

Figure 3 displays the scheme’s robustness to perturbation in initial condition. The solution 
integrated with the optimized scheme far better approximates the exact solution than the one 
employing the spectral-like approximation, even for initial data different from the one it was 
designed to resolve. This holds for both smoother and more oscillatory initial data. Although 
those examples do not give a quantitative view on the relative efficiency of the schemes for those 
initial data, one can see in both figures that by the time the solution with the optimized scheme 
developed significant error the error in the one corresponding to the spectral-like scheme is so 
large it no longer approximates the exact solution. 
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Figure 4 shows a two dimensional equation which demonstrates the robustness^" the proposed 
schemes. In this example the initial data was taken to be the Gaussian e rotated at an 

angle of Then the program searched for initial data of the form e”( niW i +n2U ^), for the integers 
1 < 7 Zi ^ 7 z 2 < 7 which yielded the best approximation to the initial data. The pentadiagonal 
schemes optimized for initial data e~ nxUJ and e~ n2W were then used to compute u x and 


respectively. In this example n\ = 3 and n 2 = 2. The resulting semi discrete system was solved 
by exact integration in Fourier space on a 32x32 grid. The plot shows a cut through the solution in 
the x direction containing the maximum point of the solution. While the solution corresponding 
to the optimized disretization closely approximates the exact solution, the solution discretized 
with the spectral-like scheme bears very little resemblance to the exact solution. 


6*1.2 Approximation of the second derivative 

The coefficients of compact schemes for various initial conditions of the form e~ au; can be found 
in Appendix B. 

Figure 5 plots absolute value of the symbol of the second derivative and the weighted error, 
for a = 2. The parameters of the optimized schemes can be found in Appendix B. The coefficients 
of the other schemes were taken from [7]. Scheme (a) is given by : 


a — yp/3 = 0, a = -jy,6 = yy,c = 0 

The coefficients of scheme (c) are : 

9 „ 696 - 1191a, 2454a - 294 1179a -344 

a = -,/3 = °,a= — ,b. — ,c 

The coefficients of the spectral-like pentadiagonal scheme (e) are : 

a = 0.50209266,/? = 0.05669169, a = 0.21564935,6 = 1.723322, c = 0.1765973 


(6.14) 


(6.15) 


(6.16) 


It can be seen that the error in the non optimized schemes is significantly larger than in the 
optimized ones. It is interesting to note that, again, for this specific data the L 2 error norm of 
the spectral-like scheme is about an order of magnitude larger than the non optimized tridiagonal 
scheme. This phenomenon suggests that the resolution efficiency is a poor estimate for discretiza- 
tions evaluation. Computing the error norms reveals that the error in the optimized tridiagonal 
scheme is about seven times smaller than in the non optimized scheme, whereas the error in the 
optimized pentadiagonal scheme is seventy times smaller than the spectral-like scheme, for this 
given data. 

The efficiency of the pentadiagonal schemes was compared by integrating the wave equation : 


d 2 U d 2 U 
dt 2 “ dx 2 


This equation was put in a system form : 



(6.17) 


(6.18) 
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This system was exact integrated on a 32 points grid and the results are given in Figures 6-7 
demonstrating the improved accuracy of the optimized scheme and its robustness, respectively. 
Figure 6 demonstrates that the optimized scheme can be used for integration time 37 times, or 
more, longer than the spectral-like scheme. In Figure 7 the scheme robustness is clearly shown for 
initial data smoother or more oscillatory than the data for which the scheme was designed. In both 
cases, by the time a significant error occurs in the solution discretized with the optimized scheme 
the solution corresponding to the spectral like scheme totally differs from the exact solution. 

The initial solution and its approximation, for the two dimensional problem in Figure 8 were 
obtained similarly to those of the example in Figure 4. While the solution integrated with the 
optimized scheme closely approximates the exact solution, it is hard to see that the solution 
corresponding to the spectral-like scheme indeed approximates the same problem. 


6*1.3 Mid cell approximation of the first derivative 


Appendix C lists the coefficients of schemes designed for various initial data. The coefficients of 
the schemes taken from [7] are listed below. Scheme (a) is given by : 

a = fi = 0, o = |(3 - 2a), b = i(22a - 1), c = 0 (6.19) 


The coefficients of scheme (c) are : 

75 37950 - 39725a 65115a -3350 25669a -6114 

a “ 354 ,/?_ °’ a_ 31368 ’ ~ 20912 62736 

The coefficients of the tenth order pentadiagonal scheme (f) are : 

96850 9675 683425 505175 69049 

° “ 288529 ’ 13 ~ 577058’° “ 865587’ 577058’° ” 11731174 


( 6 . 20 ) 


( 6 . 21 ) 


The standard compact schemes give very good resolution in this form (see Figure 9), thus, 
the improvement introduced by the optimized schemes is smaller. Optimizing the tridiagonal 
scheme yields a 6.5 smaller error norm while optimizing the pentadiagonal scheme yields a 2.5 
times smaller norm. In this case, the error norm of the optimized tridiagonal scheme is very close 
to that of the non optimized pentadiagonal scheme. 

An interesting option suggested by this approach was to optimize the ^ operator in order to 
get the best approximation for for given initial values. This has been done for the tridiagonal 
scheme which was used to integrate equation (6.17). It was compared, in Figure 10, to the 
tridiagonal scheme from [7] where both are used to approximate the second derivative in solving 
the one dimensional wave equation in the system form (6.18). Again, the optimized scheme gives 
significantly better approximation. 


6.2 Approximate time integration 

The constrained minimization (4.3)-(4.5) was solved by requiring that the solution will touch 
the stability constraint at one point while maintaining global stability and and minimizing the 
functional. The point which gives best result was found by exhaustive search. This straightforward 
approach yielded the local minima reported in this paper. Somewhat better integration schemes 
might be achieved by using more advanced optimization techniques [8]. 
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According to the general approach outlined in Section 5, one should choose an integration 
scheme which yields truncation error of similar magnitude in time and space. Since the stability 
region for several fifth order six stages explicit Runge-Kutta schemes intersects the imaginary 
axis only in a small neighborhood of the origin [5, 6] disabling time marching with large CFL, the 
optimized scheme was compared with the four stage fourth order Runge-Kutta. We preferred this 
five stage scheme which has an error norm about five times larger than the space discretization to 
the seventh order scheme which yields an error norm about eleven times smaller than the space 
discretization because of its lower computational cost. 

The analysis performed in Section 2 suggests that the integration operator should be opti- 
mized with respect to the spatial discrete operator employed, i.e., to minimize \\P(L h (u;h)At) - 
e Z,*(w/i)A*||^ j n following examples L h is the tridiagonal approximation for when initial 
data is e~ 2uj7 and a — 0.8. Appendix D contains the coefficients of integration schemes for various 
initial data when L h is the tridiagonal scheme optimized for the same initial data and a = 0.8. 

Figure 11 plots the real and imaginary parts of e Lh ^ At versus the four stage fourth order 
Runge-Kutta and the improved scheme. The norm of the imaginary part of the error was reduced 
by a factor of 31 while its real part was reduced by merely a factor of 2.3. 

Figures 12-13 shows the integration of the advection equation with those scheme on a 32 points 
grid, demonstrating the superior efficiency and robustness of the proposed schemes. In Figure 12 
one can see that the optimized scheme can be used for at least three times longer integration time 
than the Runge-Kutta scheme applied to the tridiagonal scheme from [7]. The computed error 
norms suggests the time marching error is dominant in all examples. 

The two dimensional example in Figure 14 summarizes the approach suggested in this work. 
It compares the optimized tridiagonal scheme combined with the appropriate integration formula, 
to fourth order Runge Kutta applied to non-optimized tridiagonal discretization. Although the 
analysis in Section 2 applies only to constant coefficient problems, this example shows it holds, 
heuristicly, to variable coefficient equations, as well. The initial data for this problem was obtained 
in a similar manner to that in example 4. However, instead of comparing the solutions computed 
on the 32 x 32 grid to the exact solution, they are compared to the solution on a 64 x 64 grid 
which was integrated with the optimized scheme designed for the narrowest computed Gaussian 
(a = 7). The initial data for the finer grid was obtained by bilinear interpolation from the coarser 
grid. It can be seen that the optimized scheme yields significantly more accurate solution. 

7 Conclusion 

A simple and general approach for the design of finite difference approximation of derivatives 
and integration formulas was introduced. It was used to design compact finite difference schemes 
for derivatives evaluation; and the resulting schemes were compared to previously known similar 
schemes. The guiding line has been to improve the representation of the range of wavenumbers 
appearing in the physical problem being solved, taking into account their relative amplitudes. 
This lead to an L 2 measure of the approximation. The resulting schemes combined adaptivity to 
the specific initial data by the nature of their design and robustness to perturbations in the initial 
data . The improved resolution had been demonstrated for several problems. 

A similar approach was used to design improved integration schemes, taking into account the 
spatial discretization as well as the initial data. 
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The approach suggested in this paper for optimizing discrete operators can be similarly applied 
to higher, derivatives. Its applicability to more general and complex operators should be further 
investigated. 

The use of these ideas to design boundary conditions will be presented elsewhere. 
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A Coefficients of first derivative approximations for various ini- 
tial conditions 


iiQ 

schemes with /3 = c = 0 


a = 0.3793894912, a = 1.575573790, 6 = 0.1832051925 

| 

a = 0.3534620453, a = 1.566965775, 6 = 0.1399583152 


a = 0.3461890571, a = 1.5633098070, 6 = 0.1290683071 

■besm 

a = 0.3427812069, a = 1.5614141543, b = 0.124148259 

1 

a = 0.3408027739, a = 1.5602604992, b = 0.121345048 

m 

a = 0.3395099051, a = 1.5594855939, b = 0.119534216 

KSfil 

a = 0.3385987444, a = 1.5589295176, b = 0.1182679712 


Uo 

schemes with j3 — 0 

e~ Ji 

a = 0.4303030674, a = 1.5567577428, b = 0.3451622238, c = -0.0413138317 

Bl 

a = 0.3991476265, a = 1.5636386371, b = 0.2563784492, c = -0.0217218334 


a = 0.3904091387, a = 1.5638887738, b = 0.2348222711 ,c = -0.0178927675 

ISH 

a = 0.3863287472, a = 1.5637497712, b = 0.2252138483, c = -0.0163061252 


a = 0.3839604005, a = 1.5635937780, b = 0.21976694619 , c = -0.0154399233 

■mi 

a = 0.3824122042, a = 1.5634617985, b = 0.21625718276, c = -0.0148945794 

■BM 

a = 0.3813206436, a = 1.5633544597, 6 = 0.21380659696, c = -0.0145197694 



general schemes 


a = 0.5779403671, /? = 0.0890143475 
a = 1.3030269541, 6 = 0.994883769, c = 0.0359987066 

BN 

a = 0.5801818925, (3 = 0.0877284887 
a = 1.3058941939, b = 0.9975884963, c = 0.0323380724 


a = 0.5821143744, (3 = 0.0867224075 
a = 1.3086733956, b = 0.9990906893, c = 0.0299094788 


a = 0.5831688320, f3 = 0.0862000893 
a = 1.3102698137, 6 = 0.9997174262, c = 0.0287506026 


a = 0.5838221871, /? = 0.0858844217 
a = 1.3112828763, b = 1.0000513827 , c = 0.0280789585 


a = 0.58426518608, /? = 0.0856735831 
a = 1.31197935750, 6 = 1.00025665126, c = 0.02764152958 


a = 0.58458494112, (3 = 0.08552292859 
a = 1.31248665912, b = 1.00039487751, c = 0.02733420278 
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Coefficients of second derivative approximations for various 
initial conditions 


schemes with 0 = c = 0 


a = 0.2285657609, a = 1.0139538409, b = 


a = 0.2028150072, a = 1.0598135170, b = 


a = 0.1952770765, a = 1.0716695072, 6 = 


a = 0.1917151916, a = 1.0770076313, b = 


a = 0.1896428309, a = 1.0800332355, b = 


a = 0.18828772017, a = 1.0819792783, b = 


a = 0.18733255632, a = 1.0833354275, b = 


0.4431776810 


0.3458164974 


0.3188846458 


0.3064227519 


0.29925242633 


0.29459616204 


0.29132968512 



a = 0.3125176074, a = 


a = 0.2702488609, a = 


a = 0.2580699154 , a 


a = 0.2523894606, a = 
a = 0.2491062584, a = 


a = 0.2469677390, a = 


a = 0.2454642305, a = 


schemes with 


0.7701351999, b = 


0.8863525584 , b = 


: 0.9170322739, b = 


0.9308701065,6 = 
0.9387256232, 6 = 


0.9437849227, 6 = 


0.9473144209, 6 = 


0 = 0 


0.9469577413, c = 


0.7065172637, c = 


0.6425330979, c 


0.6135153110, c : 
0.5969863585, c = 


0.5863166347, c = 


0.5788609571, c: 


-0.0920577265 


-0.0523721002 


-0.0434255409 


-0.0396064963 

-0.0374994649 


-0.0361660793 


-0.0352469171 


Uo 

general schemes 

m 

a = 0.5024750577, (3 = 0.0554440666 
a = 0.2150536435, 6 = 1.7246523136, c = 0.1761322914 

m 

a = 0.5041582074, 0 = 0.0527585356 
a = 0.2120465713, b = 1.7488409942, c = 0.1529459205 

m 

a = 0.5053986368, 0 = 0.0512444502 
a = 0.2112256102, b = 1.7609579037, c = 0.1411026601 


a = 0.5061009898, 0 = 0.0504756862 
a = 0.2110263782, b = 1.7667867767, c = 0.1353401973 

jmm 

a = 0.5065435817, 0 = 0.0500170894 
a = 0.2109783634, 6 = 1.7701652358 , c = 0.1319777431 


a = 0.5068465815, 0 = 0.0497133535 
a = 0.2109761550, 6 = 1.7723629924, c = 0.1297807226 

e~ 7jl 

a = 0.5070666579 0 = 0.0494975852 
a = 0.2109890794, b = 1.7739051293, c = 0.1282342776 
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C Coefficients of mid cell approximation of the first derivative 
for various initial conditions 


Uq 

. schemes with 0 = c = 0 

ISfll 


B&l 


BUI 

a = 0.1560892225, a = 1.0076143702, b = 0.3045640747 

fSM 

a = 0.1532174394, a = 1.0099120548, 6 = 0.2965228240 

mm 

a = 0.1515399131, a = 1.0112348225, b = 0.2918450036 

on 

a = 0.1504402935, a = 1.0120939889, b = 0.2887865980 

mm 

a = 0.1496639344, a = 1.0126967653, b = 0.2866311035 


uq 

schemes with (3 = 0 

fSM 

a = 0.2803531992, a = 0.8656018611, 6 = 0.7202754832, c = -0.0251709460 


a = 0.2421691108, a = 0.9108711860, b = 0.5897758895, c = -0.0163088538 

BH 

a = 0.2311768224, a = 0.9233491904, b = 0.5531540626, c = -0.0141496081 

warn 

a = 0.2260281312, a = 0.9290969691, 6 = 0.5361564844, c = -0.0131971911 

om 

a = 0.2230456380, a = 0.9323967450, b = 0.5263571378 ,c = -0.0126626068 

m 

a = 0.2211004185, a = 0.9345368034, b = 0.5199847302, c = -0.0123206966 


a = 0.2197316282, a = 0.9360368674, b = 0.5155096846, c = -0.0120832956 


Uq 

genera/ schemes 


a = 0.3392424034, 0 = 0.0126851467 
a = 0.7880308119, b = 0.8956208871, c = 0.0202034010 

hi 

a = 0.3364203680, 0 = 0.0159838314 
a = 0.7894607720, 6 = 0.8790559502, c = 0.0362916767 


a = 0.3359766282, 0 = 0.0164557610 
a = 0.7895453413, b = 0.8768367139, c = 0.0384827231 


a = 0.3358345755, 0 = 0.0166014190 
a0. 78955615727, b = 0.87616875736, c = 0.03914707436 

e” 5w 

a = 0.33577201042, 0 = 0.01666433833 
a = 0.78955722003, b = 0.87588406207 , c = 0.03943141540 


a = 0.33573907328, 0 = 0.01669706335 
a = 0.78955658369, b = 0.87573722427, c = 0.03957846531 

m 

a = 0.33571963682, 0 = 1.67162152050 
a = 0.78955572985, b = 0.87565178238, c = 0.03966419181 
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UQ 

schemes with ft = 0, designed to approximate 

BSM 

a = 0.2949304593, a = 0.8473898079, b = 0.7718938474, c = -0.0294227367 

mm 

a = 0.2482825125, a = 0.9037600128, b = 0.6104318128, c = -0.0176268008 

mssm 

a = 0.2349387889, a = 0.9190859222, b = 0.5656763757, c = -0.0148847202 

e -4J' 

a = 0.2287385754, a = 0.9260661646, 6 = 0.5451127700, c = -0.0137017838 

asm 

a = 0.2251628777, a = 0.9300484196, b = 0.5333228985 ,c = -0.0130455627 

mm 

a = 0.2228371850, a = 0.9326209622, b = 0.5256822919, c = -0.0126288841 

mm 

a = 0.2212037269, a = 0.9344193325, b = 0.52032910793, c = -0.0123409866 


D Coefficients of time integration scheme 


UQ 

third order schemes designed for o = 0.9 having ao = 1, — 1, i 

asm 

a 3 = 0.166281 , a 4 = 0.0397196 , o 5 = 0.0076705 

mm 

a 3 = 0.166407 , a 4 = 0.0409525 , a 5 = 0.0074510 


a 3 = 0.1664488 , a 4 = 0.04111513 , a 5 = 0.00739737 

HU 

a 3 = 0.1664805 , a 4 = 0.04121264 , a 5 = 0.00736302 


a 3 = 0.1665028 , a 4 = 0.04128218 , a 5 = 0.00733301 


a 3 = 0.1665207 , a 4 = 0.04133150 , a 5 = 0.00731074 
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Wavenumber — cj 


Wavenumber — cu 


Figure 1: Symbols (left) and absolute value of error (right) for u 0 = e 2 . (a) Sixth order 

tridiagonal scheme (/3 = c = 0) (b) Second order optimized tridiagonal scheme (/? = c = 0) (c) 
Eighth order tridiagonal scheme (j3 = 0) (d) Second order optimized tridiagonal scheme {(5 = 0) 



T = 1000. 


T = 10000. 




Figure 2: Long time integration of the equation Ut = u x , uq — e 2u> , <r — 0.8. (a) Pentadiagonal 
scheme optimized for uq = e -2u/ (b) Spectral-like pentadiagonal scheme (c) Exact solution. 


20 





T « 1800. 


T * 15000. 




Figure 3: Long time integration of the equation u* = a — 0.8. Initial solution on the left 
figure was Uo — e” w ; on the right figure it was Uo = . 

(a) PentadiagonaJ scheme optimized for uo = e~ 2u/ (b) Spectral-like pentadiagonal scheme (c) 
Exact solution. 
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T = 6500. 



Figure 4: Integration of the equation Ui — u x + Uy, <r — 0.8 using pentadiagonal schemes. Initial 
solution was u Q = rotated at an angle of \ . This data was approximated by unrotated 

gaussian € -( 3a, ?+ 2u; 2 ). (a) Optimized pentadiagonal scheme (b) Spectral-like pentadiagonal scheme 
(c) Exact solution. 
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Wavenumber — gj 


Wavenumber — cj 


Figure 5: Symbol (left) and absolute value of error (right) for ^ = e" 2u;2 . (a) Sixth order 

tridiagonal scheme (/? = c = 0) (b) Second order optimized tridiagonal scheme (j3 — c = 0) (c) 
Eighth order tridiagonal scheme (/3 = 0) (d) Second order optimized tridiagonal scheme (/? = 0) 
(e) Spectral-like pentadiagonal (f) Optimized pentadiagonal. (g) Exact symbol. Schemes were 
optimized for uq — e~ lu} . 


T * 200. 


T « 7500. 




Figure 6: Long time integration for u ti = u xx , u 0 = e 2a>2 , a = 0.8. (a) Pentadiagonal scheme 
optimized for Uo = e ~ 2uj 2 (b) Spectral-bke pentadiagonal scheme (c) Exact solution. 
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T = 3200. 


T = 16000. 




Figure 7: Long time integration for u tt = u xx , a = 0.8. Initial solution on the left figure was 
u 0 = e' 1 " 2 ; on the right figure it was u 0 = e ~ 4 “ 2 . (a) Pentadiagonal scheme optimized for 
u 0 = e~ 2u>2 (b) Spectral-like pentadiagonal scheme (c) Exact solution. 


24 





T = 6500 



0 2 4 6 

X 


Figure 8: Long time integration for u tt = u xx + u yy , a = 0.8 using pentadiagonal schemes. 
Initial solution on the left figure was u 0 = rotated at an angle of This data 

was approximated by unrotated gaussian e ~ . (a) Optimized pentadiagonal scheme (b) 
Spectral-like pentadiagonal (c) Exact solution. 
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Wavenumber— u 


Figure 9: Symbol for mid-cell discretizations of u 0 = e~ 2u>2 . (a) Sixth order tridiagonal 
scheme (/3 = c = 0) (b) Second order optimized tridiagonal scheme (/? = c = 0) (c) Eighth order 
tridiagonal scheme (/ 3 = 0) (d) Second order optimized tridiagonal scheme (/ 3 — 0) (e) Tenth 
order pentadiagonal (f) Optimized pentadiagonal. (g) Exact symbol. Schemes were optimized for 
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T = 4500 




X 

Figure 10: Long time integration for u u = u xx , uq = e ~ 2w 2 a = 0.8. (a) Non optimized tridiagonal 
scheme (b) Tridiagonal mid-cell discretization scheme of ^ optimized to approximate when 
uo = e ~ 2w 2 (c) Exact solution. 
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Figure 11: Real and imaginary parts of approximations to e L where L h (u>) is the symbol 
of the tridiagonal scheme for jj- optimized for uq — t~ 2u> and o — 0.8. (a) Five stage scheme 
optimized for the same er (b) Fourth order Runge-Kutta (c) Exact time integration. 


T - 150. 


T - 450. 




Figure 12: Integration of u t = u x , u 0 = e~ 2jl , a = 0.8. The space derivative is computed using 
the tridiagonal compact scheme optimized for the same initial date and a. (a) Five stage scheme 
optimized for this scheme and CFL (b) Fourth order Runge Kutta (c) Exact time integration 
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T = 45. 


T = 1000. 




Figure 13: Integration of u t = u x , a = 0.8. Left: tio = e ~^ ; Right: uo = e“ 4u/2 . (a) Five 
stage scheme optimized for this scheme and CFL (b) Fourth order Runge Kutta (c) Exact time 
integration. 
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T = 400. 



Figure 14: Integration of u t = u x + 0.5(1 + 0.6stn(27rj/))«j„ a = 0.8 using tridiagonal schemes. 
Initial solution was u 0 = rotated at an angle of This data was approximated 

by unrotated gaussian e~ (a) Optimized tridiagonal scheme and optimized marching 
scheme (b) Tridiagonal scheme integrated with fourth order Runge-Kutta. (c) A fine grid solution 
(practically exact) 
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